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Abstract 

Using the unintegrated gluon distribution obtained from numerical simulations of the 
Balitsky-Kovchegov equation with running coupling, we obtain a very good description of 
RHIC data on single inclusive hadron production at forward rapidities in both p+p and 
d+Au collisions. No K-factors are needed for charged hadrons, whereas for pion production 
a rapidity independent K-factor of order 1/3 is needed. Extrapolating to LHC energies, 
we calculate nuclear modification factors for light hadrons in p+Pb collision, as well as the 
contribution of initial state effects to the suppression of the nuclear modification factor in 
Pb+Pb collisions. 

1 Introduction 

The suppression of particle production at forward rapidities in d+Au collisions compared to p+p 
collisions, experimentally observed at RHIC i2j, constitutes one of the most compelling in- 
dications for the presence of non-linear QCD evolution effects in presently available data. The 
appropriate framework to study the nuclear wave function in this non-linear QCD saturation 
regime is the Color Glass Condensate (CGC), see e.g. the reviews |3l H] and references therein. 
The CGC is endowed with a set of non-linear pQCD evolution equations, the JIMWLK equa- 
tions, which in the large- A^c limit reduce to the Balitsky-Kovchegov (BK) equation [5], [6]. The 
BK- JIMWLK equations can be interpreted as a renormalization group equation for the Bjorken-x 
evolution of the unintegrated gluon distribution, and more generally of n-point correlators av- 
eraged over the nuclear wave function, in which both linear radiative processes and non-linear 
recombination effects are included. 
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Indeed, the observed reduction of the forward hadron yield in d+Au colhsions was predicted, 
albeit at a qualitative level, in [HIH], where it was directly related to the shadowing hnilt in the wave 
function of the gold nucleus due to the enhanced role of non-linear effects in its evolution towards 
larger rapidities (smaller x). Later on, a better quantitative description of the d+Au forward 
hadron yields was achieved in the CGC calculations of [9l |10l [111 [12] • These works relied on the 
use of models for the unintegrated gluon distribution of the gold nucleus inspired by approximate 
solutions of the BK equation. Relevant dynamical features in these models where either taken 
from analyses of lepton-proton scattering data or directly fitted to data, such as the anomalous 
dimension or the rapidity dependence of the saturation momentum Qg, the scale below which 
non-linear effects become important. More detailed analytical and phenomenological analyses of 
the corresponding nuclear modification factors were carried out in [13] and [HI [15] respectively. 

The reason why the BK-JIMWLK equations, the most robust theoretical tool available to 
describe the small- a; dependence of the nuclear wave function, have not been directly used in 
phenomenological studies, is that they were originally derived at leading-logarithmic accuracy 
only. It was quickly understood that this was not good enough: in analytical |T6l [T71 [TH [19] and 
numerical [201 l22l [23] studies of the original leading-order (LO) BK equation, the growth of 
the saturation scale was determined to be ~ a;~^^°, with Xj^o ~ 4.8 ^cus, incompatible with 
the phenomenology of lepton-hadron scattering which demands ~ Moreover there 

were hints that higher-order corrections would restore the compatibility of these values 

However, such insufficiency of the theory has been (at least partially) fixed through the re- 
cent calculation of the next-to-leading order (NLO) evolution equation. First, running-coupling 
corrections to the LO BK-JIMWLK kernel were derived in [25l [23 [27] . Then, the full NLO BK 
equation was obtained [28]. As demonstrated in [29], one can account for most of the higher-order 
contribution with running-coupling corrections only. Importantly, higher order corrections bring 
the evolution speed. A, down to values compatible with experimental data, among other interesting 
dynamical effects, thus narrowing the gap between theory and data. 

Indeed, first steps in promoting the BK equation with running-coupling corrections (referred 
to as rcBK henceforth) to an operational phenomenological tool have been taken in Refs. [5CT| \3T\ 
[321 [S3]- In [SD], a good description of the rapidity and collision energy dependence of the hadron 
multiplicities in Au+Au collisions measured at RHIC was achieved. Ref. [31] demonstrated the 
ability of the rcBK equation to account for the small-x behavior of the total (-F2) and longitudinal 
{Fl) structure functions measured in e+p scattering experiments. Then it was shown in [32j 
that the proton scattering amplitude fitted to data in [31] allows a good simultaneous description 
of both the proton diffractive structure function measured at HERA and of the forward hadron 
yields measured in p+p collisions at RHIC. Finally, in [SB] a good description of the few nuclear 
structure functions known at small x from e+A experiments was obtained. 

Together, these works yield a consistent picture that present experiments can probe the non- 
linear part of the hadronic and nuclear wave functions at small x, and that they can be successfully 
described by the CGC effective theory of QCD at high energies. In this work we provide a good 
description of the single inclusive hadron (charged hadron and neutral pions) yields measured 
in p+p and d+Au collisions at RHIC at forward rapidities {yh > 2), with unintegrated gluon 
distributions obtained from the rcBK equation. We also extrapolate our results to LHC energies 
and predict forward particle production in p+p, p+Pb collisions, that we present through nuclear 
modification factors for light hadrons. In the case of Pb+Pb collisions, we are able to give the 
contribution of initial state effects to the suppression of the nuclear modification factor. 
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2 Inclusive hadron spectra in d+Au collisions at RHIC 



According to Ref. [M], the differential cross section for forward hadron production in proton- 
nucleus collisions is given by 

+ Xifg/p{Xi,p^) Na(^X2,j^ Dh/g{z,p'^^) , (1) 

where pt and r/h are the transverse momentum and rapidity of the produced hadron, and fi/p 
and Dh/i refer to the parton distribution function of the incoming proton and to the final-state 
hadron fragmentation function respectively. Here we will use the CTEQ6 NLO p.d.f's [35] and 
the DSS NLO fragmentation functions [361 EZ]- In writing Eq. ([1]) we have assumed that the 
factorization and fragmentation scales are both equal to the transverse momentum of the produced 
hadron. For light hadron production discussed here, the difference between the rapidity and 
pseudo-rapidity, rjh, of the produced hadron can be neglected, yielding the following kinematics: 
Xf = A/m^ + p2/ ^snn exp {r]h) ^ Pt/V^NN exp (y^), Xi = Xp/z and X2 = Xi exp (-2?/,,), with 
^/sNN the collision energy per nucleon. Finally, the unintegrated gluon distributions (udg's) Nf(^a) 
describe the scattering of a hard valence quark (gluon) from the projectile on the saturated small-x 
glue of the target, either a nucleus or a proton. In order to avoid contamination from large (small)- 
X effects in the target (projectile), we will restrict ourselves to the study of the forward region 
Hh ^ 2 both at RHIC and LHC energies, such that Xi ^ Xq and X2 <^ Xq, where xq is the x- 
value where the small-x evolution starts (see below). Similar to previous approaches, we allow the 
possibility of a K-factor to absorb the effect of higher order corrections. For instance there is no 
a^-order term in Eq. ([T]), we shall only implement running-coupling corrections in the X2 evolution 
of Nf{a), but in principle they also affect the cross section [38j. 

The udg's Nf{A) are given by the two-dimensional Fourier transform of the imaginary part of 
the forward dipole-target scattering amplitude in the fundamental (F) or adjoint (A) representa- 
tion, Mf{a), respectively: 

NFiA){x,k) = j d^re-*^"- [l - Ar^(A)(r, y = ln(xo/x))] , (2) 

where r is the dipole size and Y is the evolution rapidity. In turn, the small-x dynamics of the 
dipole amplitudes is given by the rcBK equation: 

dUFi^jr.Y) ^ j ^2^^ ^run(^^ ^ _ y ~^ _ y^ y ~^^ 

(3) 

For simplicity, we have omitted the subscripts F{A) in the r.h.s of Eq. (|3]). Using Balitsky's 
prescription [27], the kernel in Eq. ([3]) reads 



i^^"'^(r,ri,r2) 



1 / as{rl) A 1 / as{rl) 



(4) 



where Y2 = r — ri (throughout the paper we shall use notation v = |v| for two-dimensional 
vectors) . 
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Figure 1: Negatively charged hadron and tt^ yields in proton-proton (at pseudo-rapidities (2.2, 
3.2) and (3.3, 3.8 and 4)) and deuteron-gold (at pseudo-rapidities (2.2, 3.2) and 4) collisions at 
= 200 GeV. Data by the BRAHMS and STAR collaborations. 



Following [21], we regulate the running coupling in Eqs. and (jl]) by freezing it to a constant 
value a^"^ = 0.7 in the infrared. A detailed discussion about the different prescriptions proposed 
to define the running coupling kernel and of the numerical method to solve the rcBK equation can 
be found in [29]. The only piece of information left to fully complete all the ingredients in Eq. ([1]) 
are the initial conditions for the evolution of the dipole-nucleus(proton) amplitude. Similar to 
previous works, we take them from the McLerran-Venugopalan (MV) model 



_^p(r,Y = 0) = 1 - exp 



4 V Ar 



In — 



(5) 



where Q^q is the initial saturation scale (probed by quarks), and we take A = 0.241 GeV. Contrary 
to studies of e+p data, we have discarded initial conditions a la Golec-Biernat-Wiisthoff [10], since 
their Fourier transform result in an unphysical exponential fall-off of the ugd, and therefore of the 
hadron spectra as well, at large transverse momenta. Finally, in the large- A^^c limit which we have 
implicitly assumed in order to use the rcBK equation, the gluon dipole scattering amplitude can 
be expressed in terms of the quark amplitude as 

^^A{r,Y)=2^^p{r,Y)-^^^{r,Y) . (6) 

With this setup, we obtain a very good description of RHIC data. Fig. [1] shows the comparison 
of our results with data for the invariant yield of different hadron species in p+p and d+Au 
collisions at y^s^N = 200 GeV and rapidities i/h = 2.2 and 3.2 for negative- charge hadrons (data 
by the BRAHMS collaboration [1]) and i/h = 3.3, 3.8 and 4 for neutral pions (data by the STAR 
collaboration [2]). The only free parameters adjusted to the d+Au data are xq, the value of x 
which indicates the start of the small— x evolution, and Qso, the value of the saturation scale at 
X = Xq. For the gold nucleus we obtain a quark saturation scale Q^q = 0.4 GeV^ at Xq = 0.02. 
Values of Xq between 0.015 and 0.025 are allowed within error bands, they are used to generate 
the yellow uncertainty band in Fig. [H A few comments are in order. First, the parameters 
Qso and Xq are obtained from minimum-bias data, and therefore Q^q should be considered as an 
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Figure 2: Nuclear modification factors for 7r° production in p+Pb collisions, R^pi,, for collision 
energies ^/sNN = 8.8 (left) and 6.2 TeV (right) and for rapidities yh = 2, 4, and 6. For comparison, 
the red dashed line corresponds to the same quantity calculated in the kf-factorization scheme. 



impact-parameter averaged value, the saturation scale at the center of the nucleus is bigger. We 
remind the reader that the corresponding gluon saturation scale is larger, Q^q^'""" = 0.9 GeV^. 
Second, Qso and xo are compatible with other values extracted from e+A [33j or A+A [30] data. 
They can be compared with Q^^ = 0.2 GeV^ at xq = 0.007 obtained in the case of the proton 
(in [31], Xq = 0.01 was obtained with Q^q = 0.2 GeV^). Finally, no K-factor is needed in order 
to reproduce the charged hadron yields (i.e. K = 1), whereas a rapidity independent K-factor 
K=l/3 is needed to describe the neutral pion data. Although the precise values of the K-factors 
do not have much meaning due to the 15% normalization uncertainties of the data, we do not 
have a good understanding of the strong hadron species dependence. 

3 Nuclear modification factors in p+Pb and Pb+Pb col- 
lisions at the LHC 

It is straigthforward to use formula Eq. ([1]) to calculate forward particle production in p+p and 
p+Pb at the LHC. We shall present our LHC results in terms of the nuclear modification factor 

P, = ^ / (7) 

where N^^^l is the number of binary proton-nucleon collisions in the p+Pb collision. In our 
predictions for p+Pb collisons at the LHC we use ^Vj^qH = 3.6, which is half the number of 
collisions determined in minimum bias d+Au collisions at RHIC [I]. Thus, in order to compare 
our results with experimental data one should renormalize our curves in Figs ([2]) and ([3]) to the 
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Figure 3: Nuclear modification factors for production in p+Pb collisions, -Rpp^, for collision 
energies ^/sNN = 8.8 (left) and 6.2 TeV (right) and for rapidities yh = 2, 4, and 6. For comparison, 
the red dashed line corresponds to the same quantity calculated in the kf-factorization scheme. 



number of collisions determined experimentally at the LHC Note that computing the ratio Eq. ([7]) 
removes the sensitivity to the K factors. Our Rppb calculations are displayed for two possible LHC 
energies {y/sNN = 8.8 and 6.2 TeV) in Fig. [2] for pion production and Fig. |3]for charged hadron 
production. In both cases, one observes the expected trends that Rppt decreases with increasing 
Uh, and increases with increasing pt. Our results for RpPb indicate that a significant suppression 
~ 1/2 should be expected already at not too forward rapidities. It is highly likely that the CGC 
dynamics studied here would also lead to a similar suppression at mid-rapidity (see, e.g. |I5]). 
However, to make a clear quantitative prediction for mid-rapidity one should ensure a proper 
treatment of high-x effects in the target, which is beyond the scope of this paper. 

We also compare the y = 2 and 4 curves with predictions obtained with the k^— factorization 
formalism, in order to check the validity of that approach, and especially to test up to what 
value of y it can be used. The kj— factorization formula (see Eq. ([8]) below) is valid when the 
dominant contributions to the cross section come from small values of x, for both the projectile 
(xi ^ 1) and the target (x2 ^ 1). For instance, it only includes gluonic degrees of freedom. This 
approach is clearly insufficient at very forward rapidities or large pt, where valence quarks of the 
projectile are important {xi < 1). However, as can be seen in Figs. [2] and [3], both formalisms give 
comparable results, as the lines from k^-factorization overlap with the uncertainty bands spanned 
by the results from the hybrid formalism. This seems to identify a kinematical window where both 
approximations (Eq. ([1]) and Eq. ([8]) below supplemented with parton fragmentation) are valid. 
To some extent it is not surprising that both formalisms yield comparable nuclear modification 
factors, since the suppression is ultimately rooted in the udg's themselves. 

The advantage of the kt-factorization formalism is that it can be used at mid rapidities, when 
Xi also becomes very small, invalidating formula Eq. ([1]). This is true at the LHC where at mid- 
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rapidity Xi ~ X2 ^ 1, but not at RHIC where both values of x are generally too large at ?/ = 0. 
Indeed one should keep in mind that pt c^'^/a/s/vtv and pt c'^'^/a/s/vjv are only lower values for xi 
and X2 respectively, but that through fragmentation larger values of x actually contribute more. 
The kj-factorization formula to describe inclusive gluon production reads [51] 



dr}d?pt p1 



(fb(fq ipA{xA,q,b) ipB{xB,Pt -q,Bt-b) , 



where is the impact parameter of the collision. Gluon fragmentation is not explicitely written 
down (therefore xa = Vt^-^ / \/^nn and xb = Vt^'^ / ^/snn), but it should also be accounted for. 
The unintegrated gluon distribution in Eq. ([8]), ip, is actually simply related to the one in Eq. ([1]): 

^{x,k,b) = j dhe-'^-''VlAf{r,Y = \n{xo/x),b) = eN{x,k,b) . (9) 

A detailed discussion about the definition and physical interpretation of the different udg's dis- 
cussed here can be found in [i2] . 

We had not specified the impact parameter dependence of the ugd's before because it is not 
needed in a p+p or p+A collision. Indeed in these cases one can write / (Pb fp{b) ifB^Bt — b) ^ 
ifB^Bt) J d% ipp{b). The 6-integrated proton udg does not appear in formula Eq. ([T]), it is rather 
the standard p.d.f.s that describe the (dilute) proton content in this formalism, while in the kt- 
factorization case the b dependence of (pp{b) can be safely neglected if we are only looking at ratios 
such as Eq. ([7]). As for the colhsion impact parameter Bf dependence of the target ugd ipB, 
have been dealing with minimum-bias data, therefore as mentioned before, the ugd's obtained 
with Q'^Q = 0.4 GeV^ for the gold nucleus (and 0.2 GeV^ for the proton) should be thought of Bt 
averaged udg's, but in principle we could also look at different centrality bins. Note that formula 
Eq. (El) has only been proven to be valid for p+p and p+A collisons (or dilute-dense scattering) 

HH HH HHl HG], and there are several hints that it does not hold for A+A collisions (or dense- 
dense scattering) ^7\, even if there were no final-state effects. However, numerical results seem to 
indicate that fc^— factorization breaking may not be too important in practice [30lll8lll9], but this 
could be process dependent. With the above remarks in mind, we shall use Eq. to calculate 
the initial state effects on particle production at mid-rapidity in Pb+Pb collisions at the LHC. 

We deal with the impact parameter in the following way: we assume that it factorizes as 
ipA{x,k,b) = TA{b)(pA{x,k) where Ta(&) could be for instance the Woods-Saxon profile and is 
normalized as j d% TA{b) = A. Doing so yields an impact parameter independent Raa, as the b 
integral in Eq. ([8]) is canceled by the number of collisions: 

^ / d^q <Pa{xa, q) 'fA{xB,Pt - q) ^^^^ 
/ d^q (pp{xA, q) 'pp{xB,Pt - q) 

This is acceptable for minimum bias results, and in this case the functions ifA and (fp are simply 
related to the averaged (fA and (pp, used for minimun bias p+A and p+p collisions: A^^^jj 'Pa/^p = 

^a/'^p- We shall again use A^^^jj = 3.6 in our LHC calculations. 
The nuclear modification 

RpbPb we obtain is shown in Fig. H] at the gluon level, it corre- 
sponds to gluon production immediately after the collision, i.e at proper times r = 0+. Obviously, 
in order to compare our results with data, one should convolute them with final state effects due 
to interactions of the produced gluons with the Quark Gluon Plasma, and with hadronization 
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Figure 4: Gluon level predictions from kj factorization for Pb+Pb collisions for rapidities y = 0,4. 
Solid lines correspond to an initial gluon saturation scale Qfo"""^ — ^ GeV^, and the dashed ones 
to Qf[,"°"2 = 0.8 GeV2. 

effects as well. Although this is beyond the scope of this work, our results indicate that a sizable 
part of suppression expected for single hadron yields at the LHC jSU] may be due to purely initial 
state effects. Finally, note that the value of the saturation scale probed by gluons at xq = 0.02 is 
QIq = 0.9 GeV^, this is used in Fig. IHfor minimum-bias predictions, rather a band is generated 
using = 0.8 and 1 GeV^. To study different centrality bins in Pb+Pb collisions, first one 
would have to improve our approximation for the b dependence of the ugd's. 



4 Conclusions 

In this work we have presented a good description of the hadron yields measured at forward 
rapidities in p+p and d+Au collisions at RHIC using the hybrid formalism proposed in [34j to 
describe high-energy dilute-dense scattering. The main new ingredient in our calculation is the use 
of the BK equation including running-coupling corrections to describe the Bjorken-z dependence 
of the nuclear (proton) wave functions. With the two free parameters in our work, xq and Qso 
fitted to RHIC data, we extrapolate to LHC energies without further adjustments and predict the 
suppression of the different forward hadron yields in p+Pb collisions with respect to p+p collisions. 
Using a different formalism, kj-factorization, we estimate the contribution of initial state effects to 
the nuclear modification factor in Pb+Pb collisions at the level of gluon production at the LHC. 

While our results offer an additional indication for the presence of CGC effects in RHIC data, 
the interest now focuses mostly in calibrating the expectations for the Heavy Ion program at the 
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LHC Our predictions for the LHC rely on the most up-to-date tools within the CGC formal- 
ism (Pomeron-loop corrections have been looked at but are only relevant at asymptotically 
large energies when the running coupling is included ES]) and will be useful to confirm the 
tentative conclusions reached at RHIC and to distinguish between alternative physical scenarios, 
like those proposed in [M| [55] (a complete set of predictions stemming from different formalisms 
can be found in [50]), where the suppression of forward yields at RHIC is due to the non-eikonal 
propagation of the leading parton, resulting in energy loss in the forward region. Finally, the 
proper characterization of gluon production in the early stages (before thermalization) of Pb+Pb 
collisions at the LHC would serve as crucial input for studies of the medium produced in such 
collisions, such as hydrodynamic simulations or studies of jet quenching. In the latter case, the 
suppression predicted here due to initial state effects would add to the final state effects due to 
the presence of a Quark Gluon Plasma. 
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